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As  fundamental  parameters  for  polarized-radiative-transfer  calculations,  the  single-scattering 
phase  matrix  of  irregularly  shaped  aerosol  particles  must  be  accurately  modeled.  In  this  study, 
a  scattered-field  finite-difference  time-domain  (FDTD)  model  and  a  scattered-field  pseudo- 
spectral  time-domain  (PSTD)  model  are  developed  for  light  scattering  by  arbitrarily  shaped 
dielectric  aerosols.  The  convolutional  perfectly  matched  layer  (CPML)  absorbing  boundary 
condition  (ABC)  is  used  to  truncate  the  computational  domain.  It  is  found  that  the  PSTD 
method  is  generally  more  accurate  than  the  FDTD  in  calculation  of  the  single-scattering 
properties  given  similar  spatial  cell  sizes.  Since  the  PSTD  can  use  a  coarser  grid  for  large 
particles,  it  can  lower  the  memory  requirement  in  the  calculation.  However,  the  Fourier 
transformations  in  the  PSTD  need  significantly  more  CPU  time  than  simple  subtractions  in  the 
FDTD,  and  the  fast  Fourier  transform  requires  a  power  of  2  elements  in  calculations,  thus 
using  the  PSTD  could  not  significantly  reduce  the  CPU  time  required  in  the  numerical 
modeling.  Furthermore,  because  the  scattered-field  FDTD/PSTD  equations  include  incident- 
wave  source  terms,  the  FDTD/PSTD  model  allows  for  the  inclusion  of  an  arbitrarily  incident 
wave  source,  including  a  plane  parallel  wave  or  a  Gaussian  beam  like  those  emitted  by  lasers 
usually  used  in  laboratory  particle  characterizations,  etc.  The  scattered-field  FDTD  and  PSTD 
light-scattering  models  can  be  used  to  calculate  single-scattering  properties  of  arbitrarily 
shaped  aerosol  particles  over  broad  size  and  wavelength  ranges. 

©  2013  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Remote-sensing  methods  to  determine  aerosol  properties, 
like  those  developed  for  the  NASA  Glory  mission  [1]  require 
accurate  modeling  of  the  polarized  solar  radiation  transferred 
through  the  atmosphere  composing  aerosols.  As  fundamental 
parameters  for  polarized-radiative-transfer  calculations, 
single-scattering  phase-matrix  elements  of  irregularly  shaped 
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aerosol  particles  must  be  accurately  modeled.  However,  the 
size  parameters  of  aerosols  under  the  broad  solar  spectrum 
generally  reside  in  the  range  where  neither  Rayleigh  or 
Rayleigh-Gans,  nor  geometric  optics  method  can  be  applied 
to  calculate  their  light  scattering  properties.  Although  for 
these  size  parameters  exact  algorithms  such  as  those  based 
on  Mie  theory  [2]  and  the  T-matrix  method  [3]  are  very 
efficient  at  calculating  light  scattering  from  specific,  ideal 
morphologies,  like  spheres,  spheroids,  and  cylinders,  these 
methods  generally  have  limited  applicability  to  the  real- 
world  irregularly  shaped  aerosol  particles.  Thus,  with  the 
advance  of  computing  resources,  numerical  light-scattering 
solutions  such  as  the  discrete  dipole  approximation  (DDA) 
[4-8],  the  finite-difference  time  domain  (FDTD)  technique 
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[9-16]  and  the  pseudo-spectral  time-domain  (PSTD)  method 
[17-23]  are  used  more  and  more  frequently  in  light¬ 
scattering  studies. 

The  FDTD  technique  is  formulated  by  replacing  temporal 
and  spatial  derivatives  in  Maxwell's  equations  with  their 
finite-difference  equivalences.  To  ensure  numerical  stability 
and  accurate  calculation,  the  finite-difference  spatial  cell  sizes 
are  generally  set  to  be  2d/ 20  to  2d/ 10,  where  Ad  is  the 
wavelength  inside  the  particle.  Theoretically,  the  FDTD  is  an 
accurate  method  for  any  particle  size  parameters.  However, 
for  large  particles,  the  FDTD's  requirement  for  a  small  spatial 
cell  size  results  in  formidably  large  computational  memoiy 
and  CPU  time  requirements.  Therefore,  for  large  electromag¬ 
netic  structure  scattering  problems,  the  PSTD  algorithm 
which  could  be  numerically  stable  even  with  a  spatial  cell 
size  as  large  as  2d/2  have  significant  advantages  [17],  at  least 
in  computational  memory  requirements. 

In  this  study,  we  develop  a  scattered-field  FDTD  light 
scattering  model  and  a  scattered-field  PSTD  light-scattering 
model  with  the  convolutional  perfectly  matched  layer 
(CPML)  [15,24]  absorbing  boundary  condition  (ABC)  to 
truncate  the  computational  domain,  as  a  supplementary 
development  to  the  polarized-radiative-transfer  model  for 
applications  to  irregular  aerosol  particles.  Because  the 
scattered-field  FDTD/PSTD  equations  include  incident  wave 
source  terms,  the  FDTD/PSTD  model  allows  for  the  inclusion 
of  an  arbitrarily  incident  wave  source,  including  a  plane 
parallel  wave  or  a  Gaussian  beam  like  those  emitted  by  lasers 
that  are  often  used  in  laboratory  particle  characterizations, 
etc.  In  Section  2,  the  FDTD  and  PSTD  light  scattering  models 
with  the  convolutional  perfectly  matched  layer  (CPML) 
absorbing  boundary  conditions  (ABC)  are  introduced.  In 
Section  3,  sample  numerical  results  are  shown.  Summary 
and  conclusions  are  given  in  Section  4. 


2.  Method 


2.1.  Scattered-field  finite-difference  time  domain  method 

Following  our  previous  study  for  a  light  beam's  inter¬ 
action  with  an  arbitrary  dielectric  surface  [16],  we  applied 
the  scattered-field  FDTD  technique  to  calculate  light  scat¬ 
tering  by  particles  of  arbitrary  shapes  in  free  space.  The 
scattered-field  FDTD  method  with  wave  source  terms  in  its 
update  equations  allows  more  flexibility  in  the  specifica¬ 
tion  of  the  form  of  the  incident  fields,  which  can  be  either 
plane  parallel  wave  or  a  Gaussian  beam,  etc.  [16],  Follow¬ 
ing  [16],  in  a  Cartesian  grid  system  the  x  components  of 
the  scattered  magnetic  and  electric  fields  in  the  scattered- 
field  FDTD  algorithm,  e.g.,  are  in  the  forms 

Hsxn+'/2(i,j  +  1/2,  Ac  -(- 1/2)  =  Hsxn-y2(iJ  +  1/2,  k  +  1/2) 

+  ^  x  [Esyn(lj  +  1  /2,  k  +  1  )-Esyn(i,j  +  1  /2,  k) 

+Eszn(i,j,k+l/2)-Eszn(U+  l,k+  1/2)],  (la) 


E r+1(>  +  1  /2  J,  k)  =  ( E?(i  +  1 /2  J,  k) 
j)  [H*'n+1/2(i  +  1/2, j, k-1/2) 


7  2 At /As  \  [ 
\2e  +  oAt ) 


~Hyn+,/2(i  +  1/2  j,  k  +  1/2)  +  H/n+l/2(i  +  1/2,/  +  1/2,  k) 

_Hs,n+i/2(i  +  i /2J—1  /2,  k)J 

'  [£*"+1(i  +  1/2J’ k)  +  £*"(i  +  1/2J’ k)] 

-  [£*n+1(i + 1/2J-  + y2J * k)]  -  (lb) 

where  p0  and  e0  are  the  permeability  and  permittivity  of 
free  space,  respectively;  e  and  <r  are  the  permittivity  and 
conductivity  of  the  medium  including  the  scattering  par¬ 
ticle,  respectively;  As  and  At  denote  the  spatial  cell  size 
and  time  increment,  respectively.  To  guarantee  the  numer¬ 
ical  stability  of  the  FDTD  scheme,  we  use  At=As/(2c), 
where  c  is  the  light  speed  in  free  space.  The  indices  ( i,  j,  k ) 
denote  the  center  positions  of  the  spatial  cells  in  the  FDTD 
grid.  The  time  step  is  denoted  by  integer  n.  The  positions  of 
the  magnetic  and  electric  field  components  on  a  spatial 
cell  are  identical  to  those  illustrated  in  Sun  et  al.  [13],  In 
this  study,  the  incidence  field  E'  is  set  as  a  continuous  wave 
and  analytically  given  at  any  grid  points  where  there  is 
dielectric  material. 

The  convolutional  perfectly  matched  layer  (CPML)  ABC 
developed  by  Roden  and  Gedney  [24]  is  used  to  truncate 
the  computational  domain  in  the  FDTD/PSTD  calculations. 
The  CPML  is  based  on  the  stretched-coordinate  form  of  the 
perfectly  matched  layer  (PML)  [25]  and  is  more  accurate 
than  the  uniaxial  perfectly  matched  layer  (UPML)  ABC 
[26],  The  CPML  ABC  is  more  efficient  and  more  suitable  for 
truncation  of  computational  domains  with  generalized 
materials.  For  both  scattered/total  field  formulation  and 
purely  scattered-field  formulation  of  the  FDTD/PSTD,  the 
CPML  ABC  can  be  directly  applied  to  truncate  the  calcu¬ 
lated  fields  even  where  the  incident-wave  source  exists. 
For  example,  to  match  a  CPML  along  a  boundary  to  a  lossy 
isotropic  half-space  characterized  by  permittivity  e  and 
conductivity  a  in  which  the  fields  update  equations  are 
given  as  Eqs.  (la)  and  (lb),  the  update  equations  of  Hx  and 
Ex  in  the  CPML  can  be  written  in  the  form  [15]: 

Hs.n+1/2(ij-  +  -1/2  k  +  1/2)  =  +  !  /2,  k  +  1  /2) 

-  —  { [Ezn(i,j  +  1 ,  k  +  1  /2)—Ezn(iJ,  k  +  1  /2)]  / [Ky(j  +  1  /2)As] 
Fo 

-  [Esyn(i,j  +  1/2,  k  +  1  )-Es/(iJ  +  1  /2,  k)]  /  [Kz(k  +  1  /2)As] } 

-  —  to,j  +  1/2,  k  +  1/2  )-Hsxnz(i,j  +  1/2,  k  +  1/2)1 , 

Fo  L  J 

(2) 

where 

H%(Uj  +  1/2,  k  +  1/2)  =  by(j  +  1/2  )Hsxny~\i,j  +  1/2,  k 
+1/2)  +  Cy(j  +  1  /2)[E/n(i,j  +  1,  k  +  1/2) 

- Eszn(iJ ,  k  +  1  /2)]/As,  (3a) 

+  l/2,k+l/2)  =  bz(k  +  T/2)Hs£~\i,j+  1/2,  k  +  1/2) 
+cz(k+  1  /2)[Eyn(i,j  +  1  /2,  k  +  1) 

-Eyn(i,j  +  1/2,  fc)]/As. 

Ef+1  O’  +  1  /2  J,  k)  =  ( Esxn 0  +  1  /2J,  k) 


(3b) 
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-Hz’n+1/2(i  +  1  /2,j—  1  /2,  fc)]  /  [*y(j)  As] 
-[Hf+1/2(i+l/2,j,  k+1/2) 

— Wy'n+1/2(i  +  1 12 ,j,  fc- 1  /2)]  /  [*z(k)As]  } 

-£x"+1/2(i  +  1/2,  j,  k)],  (4) 

where 

Exy+'/2(i  +  V2,  j,  k)  =  by(j)Es»-V2(i  +  1/2  J,  fc) 

+cyO)[Hf-I/2(!+  1/2J+  l/2,fc) 

— Hz’n_1/2(i  +  1  /2J-1/2,  k)]/As,  (5a) 

£^+1/2(i  +  1  /2J,  fc)  =  bz(k)E£-'/2(i  +  1  /2  J,  k) 

+cz(k)[Hyn_1/,2(i  +  1/2,  j,  fc  +  1/2) 
_Hs,n-i/2(i  +  1/2,  j,  fc— l/2)]/As.  (5b) 

The  fay,  i>z,  cy,  cz  in  these  equations  are  given  in  the  form 

(u=x,  y,  or  z): 


_  £-4t((ffu/£0KU)-Kau/£ o)) 

(6a) 

 (btl-\)au 

(6b) 

auKu  +  Kfrau 

The  CPML  properties  (ax,xx,crx) 

,  ((ly,Ky,(Jy ),  cUlCl  (CZ2,/C2,C7  %)  310 

scaled  tensor  parameters  independent  of  the  medium 
permittivity  t-  and  conductivity  c,  and  are  assigned  to  the 
FDTD  grids  in  the  CPML  in  the  form  as  [15] 


=  (x/ d )  Qx,max> 

(7a) 

Kx(X)  =  1  +  (x/d)m(/cx,max  — 

(7b) 

&x(X)  =  (X/d)m(Tx,  max- 

(7c) 

where  x  is  the  depth  in  the  CPML  and  d  is  the  CPML 
thickness  in  this  direction.  The  parameters  axmax,  xx,max 
and  (Txmax  denote  the  maximum  ax,  kx  and  ax  at  the 
outmost  layer  of  the  CPML;  e.g.,  considering  an  x-directed 
plane  wave  impinging  at  angle  9  upon  a  PEC-backed  CPML 
with  polynomial  grading  material  properties,  the  reflec¬ 
tion  factor  can  be  derived  as  [25] 


Fig.  1.  Nonzero  light-scattering  phase  matrix  elements  from  the  Lorenz-Mie  theory  (black),  the  FDTD  (red),  the  PSTD  (blue),  and  the  PSTD  with  truncation 
of  wavenumbers  (olive)  for  a  spherical  aerosol  particle  with  a  size  parameter  ofx=6  and  a  refractive  index  of  m= 1.53.  In  the  PSTD  and  FDTD  calculations, 
the  spatial  cell  size  is  set  as  1/8  incident  wavelength.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web 
version  of  this  article.) 
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R(6)  =  exp 


2 cost)  fd  , 

r 

- /  cr(x)dx 

=  exp 

eo  c  Jo 

L 

2flx,max^  60S  0 
e0C(m+  1) 


(8) 


domain  correspondent  of  the  total  field  is  obtained  from 
the  discrete  Fourier  transform  (DFT)  of  the  time  series  of 
the  field. 


Therefore,  with  a  reflection  factor  R(0)  for  normal  inci¬ 
dence,  ffx.max  can  be  defined  as 

(m  +  l)e0cln  [R(0)] 

Ox, max  —  2^j  •  w) 

As  an  accurate  approach,  R( 0)  can  range  from  10  12  to 
10-5,  /cxmax  and  axmax  can  be  chosen  with  accuracy  tests, 
and  the  numerical  errors  are  not  very  sensitive  to  them.  In 
this  study,  we  choose  R(0)=10~8,  x'Ximax=f-l,  and  ax 

max  =  0.1. 

Using  the  equations  reported  in  this  section,  a 
scattered-field  formulation  of  the  FDTD  that  includes  the 
source  term  can  be  implemented  inside  the  CPML.  The 
scattered-field  CPML  FDTD  has  great  flexibility  in  simulat¬ 
ing  electromagnetic  wave  scattering  by  dielectric  particles 
and  surfaces. 

The  single-scattering  properties  of  particles  are  calcu¬ 
lated  with  the  volume  integration  of  the  total  electric  field 
in  the  frequency  domain  [14],  The  total  electric  field  is  the 
superposition  of  the  incident  and  scattered  fields: 
E=E'+ES.  The  FDTD  simulation  is  run  for  20 Nd  time  steps, 
where  Nd  denotes  the  spatial  cell  number  in  the  largest 
dimension  of  the  computational  domain.  The  frequency 


2.2.  Scattered-field  pseudo-spectral  time  domain  method 


The  spatial  discretization  in  the  FDTD  method  causes 
dispersion  errors,  which  limit  the  spatial  cell  size  used  in 
numerical  calculations  to  not  larger  than  1/20  to  1/10  of 
the  wavelength  [13].  Therefore,  the  FDTD  method  requires 
a  large  number  of  spatial  cells  to  calculate  the  light 
scattering,  even  by  particles  of  intermediate  sizes.  The 
PSTD  algorithm  was  recently  developed  to  avoid  this 
problem,  which  ideally  needs  only  two  spatial  cells  per 
wavelength  to  discretize  the  computational  domain  and  is 
free  of  spatial  dispersion  errors  [17],  The  coarse  discretiza¬ 
tion  of  the  PSTD  algorithm  enables  this  method  to  be  an 
optimal  solution  for  large  particles.  The  errors  introduced 
in  the  PSTD  algorithm  are  claimed  to  be  only  the  temporal 
discretization  error  [18],  For  a  3-diemensional  (3D)  pro¬ 
blem,  the  PSTD  requires  a  stability  criterion  as  [18] 


At  < 


2  As 

\[2>nC 


(10) 


which  is  a  stability  criterion  with  smaller  temporal  dis¬ 
cretization  than  that  of  the  FDTD.  In  this  study,  we  set 
At=As/(3c)  for  the  PSTD  calculations. 


Fig.  2.  Same  as  in  Fig.  1,  but  in  the  PSTD  and  FDTD  calculations,  the  spatial  cell  size  is  set  as  1/16  incident  wavelength. 
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The  essence  of  the  PSTD  is  to  use  Fourier  integrals  to 
obtain  the  derivative  of  a  field  f[u) 


ff(U) 

dll 


d 

dU 


1 J  F[f(u)]e  lku dkj 

/oo 

ikF[f(u)]e~‘kudk  =  -F^1  {ifeF[f(u)]}, 

-OO 


(11) 


where  u  denotes  a  spatial  coordinate  (x,  y,  or  z);  F  and  F_1 
denote  the  forward  and  inverse  Fourier  transforms;  k 
denotes  the  wavenumber  variable  in  the  Fourier  trans¬ 
form;  and  i  =  V— 1.  This  is  a  well-known  formula  in 
mathematics  with  a  long  history  [27,28], 

Based  on  the  scattered-field  CPML  FDTD  method  intro¬ 
duced  in  Section  2.1,  we  can  obtain  the  scattered-field 
CPML  PSTD  algorithm  by  replacing  all  the  finite-difference 
approximations  of  spatial  derivatives  of  fields  in  the  FDTD 
update  equations,  including  those  in  the  CPML,  with  the 
precise  spatial  derivatives  from  the  Fourier  transforms  in 
Eq.  (11).  Since  all  the  field  components  are  calculated  at 
spatial  cell  centers  in  the  PSTD  algorithm,  the  spatial 
position  shifts  of  “+1”,  “+1/2”,  and  “-1/2”  in  all  of  the 
FDTD  update  equations,  including  those  in  the  CPML  in 
Section  2.1,  are  removed  for  the  PSTD  field  update  equa¬ 
tions.  No  other  modifications  are  necessary  in  this  FDTD- 
to-PSTD  modification.  The  PSTD  is  run  for  30 Nd  time  steps 
in  this  study. 


Note  here  that  the  forward  and  inverse  Fourier  trans¬ 
forms  are  performed  using  the  fast  Fourier  transform  (FFT) 
and  the  inverse  fast  Fourier  transform  (1FFT)  codes  given  in 
Numerical  Recipes  [29],  The  FFT  and  IFFT  are  performed 
spatially  for  each  time  step,  along  each  direction  of  x,  y, 
and  z,  using  the  field  throughout  the  whole  computational 
domain,  including  the  CPML,  to  obtain  the  spatial  deriva¬ 
tives  of  both  electric  and  magnetic  field  components  at 
each  grid  point,  prior  to  the  field  update  calculations. 

Our  practice  shows  that  the  PSTD  does  avoid  the  spatial 
dispersion  errors  that  can  cause  the  FDTD  method  to 
become  numerically  unstable  when  too  few  spatial  cells 
are  used.  However,  the  PSTD  method  does  have  significant 
numerical  errors  from  different  sources,  in  addition  to  the 
errors  caused  by  temporal  discretizations  as  claimed  by  Liu 
et  al.  [18].  The  most  pronounced  error  source  is  the  so- 
called  Gibbs  phenomenon  [18],  which  involves  both  the 
fact  that  Fourier  sums  overshoot  at  a  discontinuity,  and 
that  this  overshoot  does  not  disappear  as  the  frequency 
increases.  For  electromagnetic  wave  propagation  in  homo¬ 
geneous  media,  the  PSTD  works  ideally  because  there  is  no 
Gibbs  effect.  But  for  light-scattering  problems,  a  particular 
medium-and-field-discontinuity  case,  the  Gibbs  effect 
causes  significant  errors.  Since  the  Fourier  transforms  in 
the  PSTD  algorithm  are  conducted  throughout  the  whole 
computational  domain  at  each  time  step,  and  the  field 


0  30  60  90  120  150  180  0  30  60  90  120  150  180 

Scattering  Angle  (deg)  Scattering  Angle  (deg) 


Fig.  3.  Same  as  in  Fig.  1,  but  in  the  PSTD  and  FDTD  calculations,  the  spatial  cell  size  is  set  as  1/24  incident  wavelength. 
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derivatives  at  each  grid  point  are  actually  derived  from 
fields  throughout  the  whole  computational  domain, 
including  those  in  the  CPML,  the  Fourier  transformation 
errors  (Gibbs  effect)  at  any  field  discontinuity  are  trans¬ 
mitted  everywhere  and  to  every  time  step.  These  errors 
contaminate  the  whole  computational  domain  and  even 
cause  numerical  instability  of  the  PSTD  method  [23],  Liu 
et  al.  [18]  state  that  using  purely  scattered  fields  in  the 
calculation  could  decrease  the  Gibbs  effect,  and  this  is  also 
the  reason  why  we  use  scattered-field  FDTD  and  PSTD 
algorithm  in  this  study. 

3.  Results 

Light  scattering  is  a  typical  electromagnetic  problem 
with  discontinuities.  Though  the  PSTD  method  is  claimed 
to  be  a  rigorous  algorithm  without  numerical  dispersion 
errors,  its  Fourier  transformation  errors  (Gibbs  effect)  due 
to  these  discontinuities  and  finite  sums  of  spectral  terms 
are  actually  very  significant.  If  there  is  no  special  treatment 
to  control  these  errors,  the  original  PSTD  is  generally 
numerically  unstable  when  particle  size  parameter  is  large 
[23],  Using  volume-averaged  dielectric  constant  [23]  could 
decrease  the  discontinuities  of  the  fields,  but  could  also 
alter  the  scattering  physics  at  the  particle  edge.  Truncation 


of  the  wavenumbers  in  the  inverse  Fourier  transformation 
to  eliminate  high  wavenumber  terms  [23]  can  stabilize  the 
PSTD  algorithm,  but  can  also  cause  errors  in  the  numerical 
results  due  to  the  unphysical  manipulation  of  the  scattered 
fields.  In  this  study,  we  will  not  average  the  dielectric 
constant  in  both  the  FDTD  and  PSTD  calculations.  To  see 
the  errors  caused  by  the  wavenumber  truncation,  we  will 
compare  the  results  from  the  original  PSTD  and  the  PSTD 
with  truncation  of  the  wavenumbers  in  the  inverse  Fourier 
transformation.  We  choose  to  truncate  10%  of  the  wave- 
numbers  at  the  high  end  in  the  inverse  Fourier  transfor¬ 
mation  for  the  comparisons  in  this  study. 

Fig.  1  shows  the  nonzero  light-scattering  phase  matrix 
elements  from  the  Lorenz-Mie  theory  (black),  the  FDTD 
(red),  the  PSTD  (blue),  and  the  PSTD  with  truncation  of 
wavenumbers  (olive)  for  a  spherical  aerosol  particle  with  a 
size  parameter  of  x=6  and  a  refractive  index  of  m=  1.53  for 
dust  aerosols.  The  reason  why  we  choose  non-absorbing 
particles  with  large  refractive  index  is  that  light  scattering 
by  these  type  of  particles  is  the  most  difficult  to  be 
accurately  approached  by  numerical  models  like  the 
DDA/FDTD/PSTD.  If  a  numerical  model  is  accurate  on  these 
particles,  it  can  perform  well  on  aerosols  under  a  wide 
incident  light  spectrum.  In  the  PSTD  and  FDTD  calculations 
in  Fig.  1,  the  spatial  cell  size  is  set  as  1/8  incident 


Fig.  4.  Nonzero  light-scattering  phase  matrix  elements  from  the  Lorenz-Mie  theory  (black),  the  FDTD  (red),  the  PSTD  (blue),  and  the  PSTD  with  truncation 
of  wavenumbers  (olive)  for  a  spherical  aerosol  particle  with  a  size  parameter  ofx=12  and  a  refractive  index  of  m=1.53.  In  the  PSTD  and  FDTD  calculations, 
the  spatial  cell  size  is  set  as  1/8  incident  wavelength.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web 
version  of  this  article.) 
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wavelength.  We  can  see  that  due  to  the  big  spatial  cell  size 
of  1  / 8  incident  wavelength,  the  FDTD  results  in  significant 
numerical  dispersion  errors  in  all  of  the  light  scattering 
phase  matrix  elements.  However,  the  PSTD  results  gener¬ 
ally  follow  those  of  the  Lorenz-Mie  curves,  though  large 
errors  do  exist  at  many  scattering  angles.  When  smaller 
spatial  cell  sizes  of  1/16  and  1/24  incident  wavelength  are 
used,  as  shown  in  Figs.  2  and  3  respectively,  the  FDTD 
dispersion  errors  are  greatly  reduced.  However,  the  PSTD 
results  are  affected  by  the  spatial  cell  size  not  as  signifi¬ 
cantly  as  in  the  FDTD,  as  shown  in  Figs.  2  and  3.  This  is 
because  the  errors  of  the  PSTD  are  majorly  due  to  Gibbs 
phenomenon  and  in  nature  different  from  the  numerical 
dispersion  errors  of  the  FDTD  algorithm.  In  principle, 
numerical  dispersion  errors  can  be  totally  removed  by 
infinitely  decreasing  the  spatial  cell  size,  but  the  Gibbs 
effect  errors  of  the  PSTD  cannot.  Also,  from  Figs.  1-3  we 
can  see  that  truncation  of  wavenumbers  in  the  inverse 
Fourier  transformation  in  the  PSTD  can  cause  significant 
errors  in  the  single-scattering  phase  matrix  elements. 

For  a  larger  spherical  particle  with  a  size  parameter  of 
12  and  a  refractive  index  of  1.53,  Fig.  4  shows  that  when 
the  spatial  cell  size  is  set  as  1  / 8  incident  wavelength,  the 
PSTD  results  for  this  larger  particle  generally  follow  those 
from  the  Lorenz-Mie  theory,  though  significant  errors 
exist  at  most  scattering  angles;  whereas,  those  from  the 
FDTD  are  quite  different  from  the  exact  results  due  to  the 


numerical  dispersion  error  of  the  FDTD  method.  When 
smaller  spatial  cell  sizes  of  1/16  and  1/24  incident  wave¬ 
length  are  used,  as  shown  in  Figs.  5  and  6  respectively,  the 
FDTD  errors  are  greatly  reduced,  but  the  PSTD  errors  are 
still  not  significantly  changed.  For  the  PSTD,  the  truncation 
of  wavenumbers  cause  more  significant  errors  when  the 
spatial  cell  size  is  bigger. 

The  light  scattering  efficiencies  (Q^)  and  asymmetry 
factors  (g)  for  the  cases  in  Figs.  1-3  and  4-6  are  listed  in 
Tables  1  and  2,  respectively.  We  can  see  that  scattering 
efficiency  and  the  asymmetry  parameter  from  the  PSTD 
are  significantly  more  accurate  than  those  from  the  FDTD. 
However,  truncation  of  the  wavenumbers  in  the  PSTD  can 
cause  large  errors  in  the  scattering  efficiency  and  asym¬ 
metry  parameter. 

In  summary,  the  scattered-field  FDTD  method  with  the 
CPML  ABC  still  needs  spatial  cell  sizes  of  ~2d/20  to 
approach  a  good  accuracy  in  calculated  single-scattering 
properties.  But  a  spatial  cell  size  of  ~2d/10  in  the  PSTD  can 
result  in  accurate  results.  So  the  PSTD  is  more  suitable  for 
calculation  of  light  scattering  by  larger  particles.  However, 
the  PSTD  has  significant  Gibbs  effect  errors  which  cannot 
be  reduced  even  by  significantly  decreasing  the  spatial  cell 
size,  especially  for  large  particles.  Also,  though  the  PSTD 
requires  much  lower  spatial  resolution  (Asas/ld/10)  than 
the  FDTD  (As«/ld/20)  for  an  accurate  calculation  of  the 
single-scattering  properties  of  particles,  its  two  Fourier 


Fig.  5.  Same  as  in  Fig.  4,  but  in  the  PSTD  and  FDTD  calculations,  the  spatial  cell  size  is  set  as  1/16  incident  wavelength. 
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Scattering  Angle  (deg) 

Fig.  6.  Same  as  in  Fig.  4,  but  in  the  PSTD  and  FDTD  calculations,  the  spatial  cell  size  is  set  as  1/24  incident  wavelength. 


Table  1 

Light  scattering  efficiencies  Qs  and  asymmetry  factors  g  of  a  spherical 
particle  with  a  size  parameter  x= 6  and  a  refractive  index  m=1.53  for  the 
cases  in  Figs.  1-3.  Also  shown  are  the  relative  errors  of  the  FDTD,  PSTD, 
and  the  PSTD  with  wavenumber  truncation  (PSTDt). 


Qs 

g 

Lorenz-Mie 

2.456 

0.581 

FDTD  (A=2/8) 

1.801  (-26.6%) 

0.407  (-29.9%) 

FDTD  (A=2/16) 

2.249  (-8.4%) 

0.543  (-7.0%) 

FDTD  (A=2/24) 

2.374  (-3.3%) 

0.561  (-3.5%) 

PSTD  (A=2/8) 

2.359  (-3.9%) 

0.590  (1.5%) 

PSTDt  (A=2/8) 

3.213  (30.8%) 

0.613  (5.5%) 

PSTD  (A=2/16) 

2.554  (4.0%) 

0.594  (2.2%) 

PSTDt  (A=2/16) 

2.243  (-8.7%) 

0.594  (2.2%) 

PSTD  (A=2/24) 

2.502  (1.9%) 

0.587  (1.0%) 

PSTDt(A=2/24) 

2.409  (-1.9%) 

0.585  (0.7%) 

transformations  make  its  CPU  time  requirement  still 
comparable  to  a  FDTD  calculation  with  finer  spatial  cells. 
For  example,  our  test  shows  that  for  a  computational 
domain  of  128  x  128  x  128  spatial  cells,  the  PSTD  algo¬ 
rithm  needs  about  10  times  CPU  time  that  the  FDTD 
requires  for  the  simulation.  This  means  that  even  if  the 
PSTD  simulation  is  done  on  a  rougher  spatial  grid  like 
64  x  64  x  64  spatial  cells,  the  CPU  time  needed  in  the  PSTD 
calculation  is  still  larger  than  that  required  by  the  FDTD  for 
the  128  x  128  x  128  spatial  cells.  Note  here  that  the  FDTD 


Table  2 

Same  as  in  Table  1,  but  for  a  spherical  particle  with  a  size  parameter 
x=12. 


4 

g 

Lorenz-Mie 

2.549 

0.651 

FDTD  (A=2/8) 

1.637  (-35.7%) 

0.643  (-1.2%) 

FDTD  (A=2/16) 

2.209  (-13.3%) 

0.695  (6.7%) 

FDTD  (A =2/24) 

2.483  (-2.5%) 

0.656  (0.7%) 

PSTD  (A =2/8) 

2.193  (-13.9%) 

0.654  (0.5%) 

PSTDt  (A=2/8) 

2.695  (5.7%) 

0.715  (9.8%) 

PSTD  (A=2/16) 

2.585  (1.4%) 

0.673  (3.4%) 

PSTDt  (A=2/16) 

2.209  (-13.3%) 

0.636  (-2.3%) 

PSTD  (A=2/24) 

2.580  (1.2%) 

0.674  (3.5%) 

PSTDt  (A=2/24) 

2.574  (1.0%) 

0.700  (7.5%) 

and  PSTD  were  systematically  compared  with  the  DDA  for 
light  scattering  simulations,  although  the  FDTD  and  PSTD 
were  not  compared  directly  in  [30,31], 

4.  Conclusions 

In  this  study,  a  scattered-field  finite-difference  time 
domain  (FDTD)  and  a  scattered-field  pseudo-spectral  time 
domain  (PSTD)  model  are  developed  for  light  scattering  by 
arbitrarily  shaped  dielectric  aerosols.  The  convolutional 
perfectly  matched  layer  (CPML)  absorbing  boundary  con¬ 
dition  (ABC)  is  used  to  truncate  the  computational  domain. 
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Because  the  scattered-field  FDTD/PSTD  equations  include 
incident  wave  source  terms,  the  FDTD/PSTD  model  allows 
for  the  inclusion  of  an  arbitrarily  incident  wave  source, 
including  a  plane  parallel  wave  or  a  Gaussian  beam  like 
those  emitted  by  lasers  usually  used  in  laboratory  particle 
characterizations,  etc.  Numerical  results  show  that  single¬ 
scattering  properties  of  spherical  particles  from  the 
scattered-field  CPML  FDTD  are  close  to  those  from 
Lorenz-Mie  theory  only  when  small  spatial  cell  sizes  are 
used.  Flowever,  the  PSTD  can  produce  more  accurate 
results  when  spatial  cell  size  is  large.  As  a  trade-off  for 
numerical  stability,  truncation  of  wavenumbers  in  the 
inverse  Fourier  transformation  in  the  PSTD  can  cause 
significant  errors  in  the  calculated  single-scattering  prop¬ 
erties  of  aerosols,  especially  in  scattering  efficiencies.  Same 
as  the  FDTD,  the  PSTD  has  no  preference  to  a  particle 
shape,  though  different  particle  shapes  can  numerically 
cause  some  difference  in  results  due  to  the  imperfect  ABC 
and  Gibbs  effect.  Examples  of  application  of  the  PSTD  to 
nonspherical  particles  can  be  found  in  [18], 

The  PSTD  can  be  a  good  algorithm  for  coarse-mode 
aerosols  that  requires  much  less  computing  memory,  though 
two  Fourier  transformations  and  a  power  of  2  elements 
required  by  the  fast  Fourier  transform  make  its  CPU  time 
requirement  still  comparable  to  a  FDTD  calculation  which 
requires  fine  spatial  cubic  cells.  Therefore,  depending  on  the 
computer's  available  memory  and  CPU,  the  FDTD  method 
could  be  applied  for  small  aerosol  particles,  and  for  larger 
particles  the  PSTD  method  can  be  used.  Employing  the 
advantages  of  both  methods,  single-scattering  properties  of 
arbitrarily  shaped  aerosol  particles  can  be  calculated  over 
broad  size  and  wavelength  spectra. 
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